//#define fast_norm3( poot ) (poot).w = native_rsqrt(mad((poot).x, (poot).x, mad((poot).y, (poot).y, (poot).z*(poot).z))); (poot).x *= (poot).w; (poot).y *= (poot).w; (poot).z *= (poot).w; (poot).w = 0.0f;
//#define fast_norm3_t( temp, poot ) (temp) = native_rsqrt(mad((poot).x, (poot).x, mad((poot).y, (poot).y, (poot).z*(poot).z))); (poot).x *= (temp); (poot).y *= (temp); (poot).z *= (temp);

#define fast_norm3( v ) (v).w = native_rsqrt((v).x*(v).x + (v).y*(v).y + (v).z+(v).z); (v).x *= (v).w; (v).y *= (v).w; (v).z *= (v).w; (poot).w = 0.0f;
#define fast_norm3_t( temp, v ) (temp) = native_rsqrt((v).x*(v).x + (v).y*(v).y + (v).z*(v).z ); (v).x *= (temp); (v).y *= (temp); (v).z *= (temp);

#define fast_dot3( v1, v2 ) ((v1).x*(v2).x + (v1).y*(v2).y + (v1).z*(v2).z)
#define fast_sub3( r, v1, v2 ) (r).x = (v1).x - (v2).x; (r).y = (v1).y - (v2).y; (r).z = (v1).z - (v2).z;
#define fast_add3( r, v1, v2 ) (r).x = (v1).x + (v2).x; (r).y = (v1).y + (v2).y; (r).z = (v1).z + (v2).z;
#define fast_min3( r, v1, v2 ) (r).x = fmin((v1).x, (v2).x); (r).y = fmin((v1).y, (v2).y); (r).z = fmin((v1).z, (v2).z);
#define fast_max3( r, v1, v2 ) (r).x = fmax((v1).x, (v2).x); (r).y = fmax((v1).y, (v2).y); (r).z = fmax((v1).z, (v2).z);

#ifndef BLOCK_X
#error BLOCK_X (local size) undefined
#endif

#ifndef BLOCK_Y
#error BLOCK_Y (local size) undefined
#endif

#define NUM_THREADS ( BLOCK_X * BLOCK_Y )

typedef struct {
	float x;
	float y;
	float z;
} float3h;

int boxIntersect( float3h min1, float3h max1, float3h v, float boxSize ) {
	
	return ( min1.x < v.x + boxSize ) && \
	( min1.y < v.y + boxSize ) && \
	( min1.z < v.z + boxSize ) && \
	( max1.x > v.x - boxSize ) && \
	( max1.y > v.y - boxSize ) && \
	( max1.z > v.z - boxSize );
	
	/*if ( min1.x > v.x + boxSize ) return 0;
	if ( min1.y > v.y + boxSize ) return 0;
	if ( min1.z > v.z + boxSize ) return 0;
	if ( max1.x < v.x - boxSize ) return 0;
	if ( max1.y < v.y - boxSize ) return 0;
	if ( max1.z < v.z - boxSize ) return 0;
	return 1;*/
}

#define MASK_THREAD ( (NUM_THREADS) - 1 )

#if NUM_THREADS == 1
	#define DIV_THREAD 0
#elif NUM_THREADS == 2
	#define DIV_THREAD 1
#elif NUM_THREADS == 4
	#define DIV_THREAD 2
#elif NUM_THREADS == 8
	#define DIV_THREAD 3
#elif NUM_THREADS == 16
	#define DIV_THREAD 4
#elif NUM_THREADS == 32
	#define DIV_THREAD 5
#elif NUM_THREADS == 64
	#define DIV_THREAD 6
#elif NUM_THREADS == 128
	#define DIV_THREAD 7
#elif NUM_THREADS == 256
	#define DIV_THREAD 8
#endif

#if PIXELS_X_PER_TILE == 8
	#define SHIFT_TILE_X 3
#elif PIXELS_X_PER_TILE == 16
	#define	SHIFT_TILE_X 4
#elif PIXELS_X_PER_TILE == 32
	#define	SHIFT_TILE_X 5
#elif PIXELS_X_PER_TILE == 64
	#define	SHIFT_TILE_X 6
#elif PIXELS_X_PER_TILE == 128
	#define	SHIFT_TILE_X 7
#elif PIXELS_X_PER_TILE == 256
	#define	SHIFT_TILE_X 8
#elif PIXELS_X_PER_TILE == 512
	#define	SHIFT_TILE_X 9
#endif

#if PIXELS_Y_PER_TILE == 8
	#define SHIFT_TILE_Y 3
#elif PIXELS_Y_PER_TILE == 16
	#define	SHIFT_TILE_Y 4
#elif PIXELS_Y_PER_TILE == 32
	#define	SHIFT_TILE_Y 5
#elif PIXELS_Y_PER_TILE == 64
	#define	SHIFT_TILE_Y 6
#elif PIXELS_Y_PER_TILE == 128
	#define	SHIFT_TILE_Y 7
#elif PIXELS_Y_PER_TILE == 256
	#define	SHIFT_TILE_Y 8
#elif PIXELS_Y_PER_TILE == 512
	#define	SHIFT_TILE_Y 9
#endif

// divide by the number of threads, add 1 if there is a remainder

__kernel void simpleGenLightBuffer(	__read_only	 image2d_t	normSpecImage,
									__read_only	 image2d_t	depthImage,
									__write_only image2d_t	diffuseImage, 
									__write_only image2d_t	specularImage,
									const	 float		specularPowerScale,
									const	 float		Q1, // 2 * zNear / ( right - left )
									const	 float		Q2, // 2 * zNear / ( top - bottom )
									const	 float		C,  // - (zFar + zNear) / (zFar - zNear)
									const	 float		D,  // - (2 zFar zNear) / (zFar - zNear)
									sampler_t	nearestSampler,
									const		int2		imageSize,
									const		float		animation_time,
									__read_only image2d_t	perBlockImage1, 
									__read_only image2d_t	perBlockImage2,
									__global	int		*lightBins,
									__global	int		*lightsPerBin,
									__global	float4	*lights,
									int			binOffset,
									int			numTilesX)
{



	const int2 coordsi = { get_global_id(0), get_global_id(1) };
	const float2 coords =  { (float)coordsi.x, (float)coordsi.y };
	
	__local uint block_culled[NUM_THREADS];
	__local float4 sharedLights[NUM_THREADS];
		
	/* per block attributes */
	__local float block_cone_r2;
	__local float3h block_min, block_max;
	__local float3h block_normal;
	__local int numLights;
	__local float3h blockPos;
	__local int binID;
	__local int binStart;
	__local float2 inv_image_size;
	
	__local int fastPath;
	/* end per block attributes */

	float debug1 = 0.0f;
	float debug2 = 0.0f;
	
	float temp;

	
	const uint linear_threadID = get_local_id(0) + (BLOCK_X) * get_local_id(1); // too big for mul24
	
	float3h diffuse	= { 0.0f, 0.0f, 0.0f };
	float3h spec	= { 0.0f, 0.0f, 0.0f };

	float d1 = -native_recip(Q1);
	float d2 = -native_recip(Q2);

#if 1
		// BLOCKWISE COMPUTATIONS

	if ( linear_threadID == 0 ){
	
		inv_image_size.x = native_recip((float)imageSize.x);
		inv_image_size.y = native_recip((float)imageSize.y);
	
		int binID	= (coordsi.x >> SHIFT_TILE_X) + mul24((coordsi.y >> SHIFT_TILE_Y), numTilesX);
		binStart  = mul24(binID , binOffset);
		numLights = lightsPerBin[binID];

		fastPath = numLights < 16;

		if (!fastPath) {
	
			float2 block_coords = { get_group_id(0) , get_group_id(1) };

			float4 temp_lookup2 = read_imagef(perBlockImage2, nearestSampler, block_coords);
			float depth1 = temp_lookup2.x;
			float depth2 = temp_lookup2.y;

			block_min.x  = 2.0f * (block_coords.x * (float)(BLOCK_X))			* inv_image_size.x - 1.0f;
			block_min.y  = 2.0f * (block_coords.y * (float)(BLOCK_Y))			* inv_image_size.y - 1.0f;
			block_max.x  = 2.0f * (block_coords.x * (float)(BLOCK_X+1) - 1.0f)	* inv_image_size.x - 1.0f;
			block_max.y  = 2.0f * (block_coords.y * (float)(BLOCK_Y+1) - 1.0f)	* inv_image_size.y - 1.0f;
			
			// at this point everything is in normalized device coordinates
						
			block_min.z = -D / ( depth2 + C );
			block_max.z = -D / ( depth1 + C );
			block_min.x *= d1 * ( block_min.x < 0.0f ? block_min.z : block_max.z );
			block_max.x *= d1 * ( block_max.x < 0.0f ? block_max.z : block_min.z );
			block_min.y *= d2 * ( block_min.y < 0.0f ? block_min.z : block_max.z );
			block_max.y *= d2 * ( block_max.y < 0.0f ? block_max.z : block_min.z );

			//Get the block normal
			float4 temp_lookup1;
			temp_lookup1 = read_imagef(perBlockImage1, nearestSampler, block_coords);
			block_normal.x = temp_lookup1.x;
			block_normal.y = temp_lookup1.y;
			block_normal.z = temp_lookup1.z;
			block_cone_r2  = temp_lookup1.w;
			
			//float3h blockPos;
			blockPos.x = (block_max.x + block_min.x) * 0.5f;
			blockPos.y = (block_max.y + block_min.y) * 0.5f;
			blockPos.z = (block_max.z + block_min.z) * 0.5f;
			float shift_x = 0.5f * (block_max.x - block_min.x);
			float shift_y = 0.5f * (block_max.y - block_min.y);
			float shift_z = 0.5f * (block_max.z - block_min.z);
			blockPos.x -= shift_x * (block_normal.x > 0.0f ? 1.0f : -1.0f);
			blockPos.y -= shift_y * (block_normal.y > 0.0f ? 1.0f : -1.0f);
			blockPos.z -= shift_z * (block_normal.z > 0.0f ? 1.0f : -1.0f);
		
		}

	}	
	barrier(CLK_LOCAL_MEM_FENCE);
	
#endif

	//Compute per pixel values that are not per light
	// compute eye-sace position
	float depth = 2.0f * read_imagef(depthImage, nearestSampler, coords).x - 1.0f;
	
	float3h worldPos;
	worldPos.z  = native_divide(-D , ( depth + C ));
	worldPos.x  = (2.0f * coords.x * inv_image_size.x - 1.0f) * worldPos.z * d1;
	worldPos.y  = (2.0f * coords.y * inv_image_size.y - 1.0f) * worldPos.z * d2;

	// compute normalized eye position
	float3h eyeNormalized = worldPos;
	fast_norm3_t(temp, eyeNormalized);

	// lookup normal / specular
	float4 normalSpec		= read_imagef(normSpecImage, nearestSampler, coords);
	float specularPower		= normalSpec.w;
	// compute normal
	float3h normal			= { normalSpec.x - 0.5f, normalSpec.y - 0.5f, normalSpec.z - 0.5f};
	fast_norm3_t(temp, normal);
	float specFactor = specularPower * specularPowerScale;
	
	#define boxSize 0.2f
	#define box_size_sq (boxSize * boxSize)

#if 1
	
	int lightIndex = linear_threadID;
		
	for (int i=0; i<numLights; i++) {

		int mask = i & MASK_THREAD;
		if ( mask == 0 ) {
		
			int isCulled = 1;

			if ( lightIndex < numLights ) {
			
				isCulled = 0;
				
			
				int index = lightBins[binStart + lightIndex];
				float4 lightRead = lights[index];
								
				sharedLights[linear_threadID] = lightRead;

				if (!fastPath) {

					float3h light = { lightRead.x, lightRead.y, lightRead.z };
					bool isPoint = lightRead.w > 0.99f;
					
					if (isPoint) {
						// BOX TEST
						isCulled = !boxIntersect( block_min, block_max, light, boxSize );
					}
					
					if (!isCulled && block_cone_r2 > 0.0f) {
						// CONE TEST
						float3h lightDirection;
						if ( isPoint) {
							fast_sub3(lightDirection, light, blockPos );
						}
						else {
							lightDirection.x = light.x;
							lightDirection.y = light.y;
							lightDirection.z = light.z;
						}
						fast_norm3_t(temp, lightDirection);
						float avg_dot		= -fast_dot3(lightDirection, block_normal);
						float avg_dot_sqr	= avg_dot * avg_dot;
						float avg_dot_cub	= avg_dot_sqr * avg_dot;
						isCulled = (avg_dot > 0.0f) && ((avg_dot_sqr - avg_dot_cub) > (0.5f * block_cone_r2));
					}
				}
			}				
			// store culling results to shared memory
			block_culled[linear_threadID] = isCulled;
			barrier(CLK_LOCAL_MEM_FENCE);
			lightIndex += NUM_THREADS;
		}
			
		// if the light is culled, continue
		// we don't need to check that it's a valid index
		// because it will be culled if it's invalid
		if ( !block_culled[mask] ) {
								
			float4 lightRead = sharedLights[mask];
			int isPoint = lightRead.w == 1.0f;

			float3h light = { lightRead.x, lightRead.y, lightRead.z };
			
			float3h lightColor;
			
			float attenuation;
			if ( isPoint ) {
				fast_sub3(light, light, worldPos );
				float d_sq		  = fast_dot3(light, light);
				float f			  = native_divide(d_sq , box_size_sq );
				attenuation		  = 1.0f * fmax( 0.0f, 0.03f / d_sq * ( 1.0f - f * f ) );
				if ( attenuation > 3.0f ) attenuation = 3.0f;
				fast_norm3_t(temp, light);
				lightColor.x = 1.0f;
				lightColor.y = 0.5f;
				lightColor.z = 0.1f;
			}
			else {
				fast_norm3_t(temp, light);
				lightColor.x = 1.0f;
				lightColor.y = 1.0f;
				lightColor.z = 1.0f;
				attenuation = 0.3f;
			}
			
			// compute diffuse component
			float kd = ( attenuation * fmax( 0.0f, fast_dot3(normal, light) ) );
			diffuse.x				+= lightColor.x * kd;
			diffuse.y				+= lightColor.y * kd;
			diffuse.z				+= lightColor.z * kd;
			
			// compute specular
			float3h halfV;
			fast_sub3( halfV, light, eyeNormalized);	
			fast_norm3_t(temp, halfV );
			float halfVdotN			= fmax(0.0f, fast_dot3(halfV, normal));
			float ks				= attenuation * native_powr(halfVdotN, specFactor );
			spec.x					+= lightColor.x * ks;
			spec.y					+= lightColor.y * ks;
			spec.z					+= lightColor.z * ks;
			
			debug1 += 1.0f / 16.0f;

		}
															
	}
#endif
	
	debug2 =  numLights / 16.0f;
	
	float4 out1 = { diffuse.x, diffuse.y, diffuse.z, debug1 };
	float4 out2 = { spec.x, spec.y, spec.z, debug2 };
	
	bool validPixel = (coordsi.x < imageSize.x && coordsi.y < imageSize.y);

	if ( validPixel ) {
		// write outputs
		write_imagef( diffuseImage,  coordsi, out1 );
		write_imagef( specularImage, coordsi, out2 );
	}

}